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Abstract 

We study numerically stabilized solutions of the two-dimensional Schrodinger equa- 
tion with a cubic nonlinearity. We discuss in detail the numerical scheme used and 
explain why choosing the right numerical strategy is very important to avoid mis- 
leading results. We show that stabilized solutions are Townes solitons, a fact which 
had only been conjectured previously. Also we make a systematic study of the pa- 
rameter regions in which these structures exist. 
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1 Introduction 



In this paper we study systems ruled by the two-dimensional nonlinear Schro- 
dinger equation (NLSE) [1,2] with a cubic time-dependent nonlinearity. More 
precisely we look for solutions of the Cauchy problem 



du 



\A + g(t)\u\ 2 



u (la) 
u(x,0) = u (x) E H\R 2 ) (lb) 



where u(x, t) : IR 2 x M + — ► C, A = d 2 /dx\ + d 2 jdx\ and g(t) is a real function 
(the nonlinear coefficient) so that if g < the nonlinearity is attractive whereas 
for g > the nonlinearity is repulsive. 

When g is a real constant Eq. (la) is the cubic NLSE, which is one of the most 
important models of mathematical physics, with applications in very different 
fields [1,2]. 
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It is well known that for g < if N = 1 1 Hi? * s above a threshold value N c , 
solutions of Eq. (la) can self-focus and become singular in a finite time. This 
phenomenon is called wave collapse or blowup of the wave amplitude. More 
precisely, there is no blowup when N < N c but for any e > 0, there exist 
solutions with N = N c + e for which there is blowup [3,4]. 

Eq. (la) admits stationary solutions of the form u(x,t) = e^Q^x), where 
^n(x) verifies 

A^-2/^-2<?|^| 2 $ M = 0. (2) 

As it is precisely stated in [1], when g is negative, for each positive \i there 
exists only one solution of Eq. (2) which is real, positive and radially symmetric 
and for which / \^ p \ 2 d 2 x has the minimum value between all of the possible 
solutions of Eq. (2). Moreover, the positivity of fi ensures that this solution 
decays exponentially at infinity. This solution is called the ground state or 
Townes soliton. We will denote this solution as R^{r) which satisfies 

Ai? M - 2/iR^ - 2gRl = (3) 
Bmi2»=0, # M (0) = 0. (4) 

Once the value of \x is fixed, the power and width of are given by iV M = 
/ \R^\ 2 d 2 x and = (J \R^\ 2 r 2 d 2 x) 1 / 2 respectively. By applying scaling trans- 
formations to R^{r) it is possible to build a family of Townes solitons having 



the same shape and norm but different widths 

R,(r) = n 1/2 Ri(» 1/2 r), W, = W ± /^ 2 . (5) 
The equation which verifies the normalized soliton R^ N (r) = R^(r)/N^ 2 is 

AR^ N - 2fiR^ N - 2gN ll Rl N = 0. (6) 



From the theory of nonlinear Schrodinger equations it is known that the 
Townes soliton has exactly the critical power for blowup iV c , therefore, it sepa- 
rates in some sense the region of collapsing and expanding solutions. Moreover, 
the Townes soliton is unstable, i.e. small perturbations of this solution lead to 
either expansion of the initial data or blowup in finite time. 

The case where g is not constant but a continuous periodic function of t has 
arisen recently in different fields of applications of Eq. (la). The intuitive 
idea is that (oscillating) bound states could be obtained by combining cycles 
of positive and negative g values so that after an expansion and contraction 
regime the solution could come back to the initial state. In this way some sort 
of pulsating trapped solution, i.e., a breather, could be generated. 

This idea was first proposed in the field of nonlinear Optics [5]. In that con- 
text, two-dimensional solitary waves have been proposed for optical systems 
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[5,6,7]. There, a spatial modulation of the Kerr coefficient (the nonlinearity) 
of the optical material is used to prevent collapse so that the beam becomes 
collapsing and expanding in alternating regions and is stabilized in average. 
The same idea has been used in the field of matter waves in Refs. [8,9]. In 
Ref. [10] some general results are provided for generic nonlinearities. Finally 
in Ref. [11] stabilized vector solitons have been proposed and studied. 

The aim of this paper is to study in more detail the stabilization mechanism. 
This paper is organized as follows: First, in Sec. 2 we present the numerical 
method used to integrate the equations, which involves a Fourier pseudospec- 
tral scheme, a split-step scheme and the use of an absorbing potential which 
allows to get rid of the radiation. In Sec. 3 we discuss the phenomenon of 
stabilization of solutions of the NLSE and show quantitatively that the struc- 
ture that remains stabilized is a Townes soliton, thus confirming the conjecture 
raised in Ref. [10]. We also study the robustness of the stabilized soliton under 
parameter variations. Finally, in Sec. 4 we summarize our conclusions. 



2 Numerical methods 

2. 1 Fourier pseudoespectral scheme 

In this paper we study Eq. (la) numerically. To this end we have devel- 
oped a Fourier pseudospectral scheme for the discretization of the spatial 
derivatives combined with a split-step scheme to compute the time evolu- 
tion. Split-step schemes are based on the observation that many problems 
may be decomposed into exactly solvable parts and on the fact that the 
full problem may be approximated as a composition of the individual prob- 
lems. For instance, the solution of partial differential equations of the type 
dtu(x,t) = N (t, x, u, d x , ....) u = (A + B)u can be approximated from the 
exact solutions of the problems d t u = Au, and d t u = Bu. 

To approximate properly the solutions of Eq. (la) we include an absorbing 
potential in the border of the simulation region to extract the outgoing radi- 
ation and provide some sort of local transparent boundary conditions. This 
potential may be included in Eq. (la) as a complex term [12,13] so that the 
new equation we will discretize is 



where V a (x) is a positive real function which must be chosen to maximize the 
absorption of the radiation while minimizing the effect on trapped structures 
(i.e. to mimic the behavior of true transparent boundary conditions). More- 



.dU 



r i 



-A + g(t)\U\ 2 -iV a \U, 



(7) 
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over, U(x,t) ~ u(x, t) will satisfy specific boundary conditions on the border 
of the simulation region Q h . As fl^ — > 1R 2 and \4 — > we expect U — > it. 

Let us decompose the evolution operator in Eq. (7) by taking [14,15] 



A = ^A, (8a) 
£ = -^(0|t/| 2 -K- (8b) 

To proceed with split-step type methods it is necessary to compute the explicit 
form of the operators e^~ to ^ A and e^~ to ^ B . To get the action of the operators 
we solve the subproblems 



d t U= l -AU, (9a) 
d t U = -ig(t)\U\ 2 U -V a U. (9b) 

Since U G H 1 (M 2 ), Eq. (9a) can be solved in Fourier space. Denoting the 
spatial Fourier transform of U by U(k,t) = J-{U{x,t)) we get 

d t U = - l -k 2 u, (10) 

whose explicit solution is U(k,t) = U(k, t )e~ lk2 ^~ to ^ 2 . 

To solve (9b) we write U = p x l 2 e 1 ^ and substituting in Eq. (9b) we get, from 
the real part, p = —2V a p whose solution is 

p(t)=e- 2V ^p(t ). (11) 

The imaginary part of Eq. (9b) is <fi = —g(t)p whose solution is <f>(t) = (fi(t ) — 
p(t ) // () g(s)e~ 2Va( - s ~ u ^ds. Finally, we can write 

17(f) = U(t ) exp (-V a (t - t ) - ip(t ) j\{s)e- 2V ^ s ^ds^ . (12) 

In this way, by defining the time step r as r = (t — 1 ) / c, c G M, after a suitable 
renaming we obtain the explicit form of the operators 



e (t-t Q )A = e crA = jr-i exp (_ ?cr A; 2 /2)^, (13a) 

£ (t-to)B = e crB = exp (^_ VaCT _ J' +CT g(s)e' 2V ^ S -^ds^j . (13b) 

Specifically, when g(t) is given by g(t) = g + gicos(fit), as we have used in 
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Fig. 1. Reflection (R) and transmission (T) coefficients (see text) for the collision 
of a Townes soliton with an absorbing potential V a = Ae~( r ~ n> } / 2u V (a) T and R 
as a function of A for w a = 0.5. (b) T and R as a function of w a for A = 5. 

our simulations it is possible to write explicitly the integral in (13b) obtaining 

rt+cr 



X 



i; 

{e- 2VaCT [Q sin Q(t + cr) - 2V a cos Q{t + cr)} - (Q sin Qt - 2V a cos ttt) } . 



, (s) e-<-.» ds = (e— - 1) - 



Having the explicit form of the solutions of the subproblems (9) allows us to 
obtain the solution of Eq. (7) to any degree of accuracy. We have used two 
different splittings: the classical method of second order and a fourth order 
optimized method [16] whose equations are 



U(x, t + r) = e TA/2 e TB e rA/2 U{x, t) + C(r 3 ), (14a) 

U^X t + t) — e a ^ rA e b ^ rB e a 2 TA e b 2 TB e a 3 TA e b 3 TB e a ^^ e hrB e a 3 rA x 

x e b2r V 2r V iri V irA l7(;r, t) + C(r 5 ). (14b) 

withai = 0.0829844064174052, 6i = 0.245298957184271, a 2 = 0.396309801498368, 
b 2 = 0.604872665711080, a 3 = -0.0390563049223486, b 3 = 1/2 - (h + b 2 ), 
0,4 = 1 — 2(ai + a 2 + 03). Both methods are symmetric compositions of op- 
erators. However, method (14a) requires only one evaluation of e A and e B 
(instead of two) per step because it is possible to concatenate terms using the 
First Same As Last (FSAL) property. For method (14b) the computational 
cost per step is much higher since even using the FSAL property the number 
of evaluations of e A and e B is six. 



These schemes have many advantages. From the practical point of view the cal- 
culation of the Fourier transform, which is the most computer time-consuming 
step in the calculations, may be done by using the fast Fourier transform 
(FFT). Thus the computational cost of the method is of order O [M 2 log Af), 
being M the points number in each spatial direction of the grid which is 
quite acceptable. The use of discrete transforms to represent the continuous 



5 



Fourier transform in (14) implicitly imposes periodic boundary conditions on 
U. However, since U is expected to be negligible on the boundaries (other- 
wise the computational domain must be enlarged or the absorbing potential 
corrected) this is not an essential point. Another convenient property of these 
schemes is their preservation of the L 2 -norm of the solutions. 

2.2 The absorbing potential 

The implementation of transparent boundary conditions is essential to avoid 
misleading results as it will be clear later. This is because of the decomposition 
of initial data into trapped states plus some outgoing radiation, a situation 
which is typical of NLS problems with fixed coefficients [17,18] and also hap- 
pens in the case at hand. 

Specifically, we have chosen as absorbing potential V a = Ae~ ( - r ~ r °^/ 2w ", where 
A, r and w a characterize the absorber. The choice of the parameter values 
must be done to maximize the absorption of the outgoing radiation. To choose 
optimal values we have made several tests by making a Townes soliton collides 
with absorbing potentials of several amplitudes and widths. We have computed 
the reflection (R) and transmission (T) coefficients defined as R = \ub\/\u \, 
T = \uf\/\u \ where u , Ub and Uf are the initial data, the part of the initial 
data which moves backwards after reflection, and the part of the initial data 
which moves forward after crossing the absorbing potential, respectively. 

Typical examples of the results are shown in Fig. 1. One must choose a value 
for A large enough to ensure a good absorption, whereas the width w a of the 
absorbing potential cannot be chosen too large, since then either the effect on 
stabilized structures could be non-negligible or the integration region should 
be enlarged to non-practical sizes. 

When the absorbing potential is absent the numerical simulations are mislead- 
ing and the trapping effects are altered, because the radiation interacts with 
the stabilized structure leading to spurious destabilization. In Figs. 2 and 3 
we illustrate this behavior starting with a Townes soliton as initial data. 



3 Stabilization of solitons 

3. 1 What is the shape of stabilized solitons ? 

In Refs. [5,6,7,8,9,10] it has been shown that if g is an adequate periodic func- 
tion g(t) it is possible to obtain a stabilized structure starting from different 
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Fig. 2. Results of numerical simulations of Eq. (la) showing evolution of the initial 
data u(x,0) = R^n^x) with fi = 0.5 and g = —0.5 (W(0) = 1.09) for parameter 
values go = —2tt, g\ = 8ir, f2 = 40, without absorbing potential, (a) Evolution of the 
width W(t) and of the amplitude max x \u(x, t)\. The insets show pseudo-color plots 
of |u(x,t)| 2 for different times, (b) Evolution of the norm lino- 
types of initial data, such as Townes solitons or Gaussian functions. In Figs. 
3 and 4 we summarize results taking as initial data a Townes soliton and a 
Gaussian function, respectively. 

In both cases we get stabilization but in the case of Gaussian initial data, a 
readjustment is produced by ejecting outside the central region a significant 
part of the wave packet. This happens during the first 10 time units where the 
width increases significantly due to the contribution of the outgoing wave. This 
wave is dissipated when hits the absorbing region leading to the step in the 
norm evolution shown in Fig. 4(b). It was also conjectured in Ref. [10] that the 
stabilized solution after this process might be a scaled Townes soliton (which 
would explain why the stabilization of the Townes soliton looks smoother). In 
this section we support this conjecture quantitatively with precise numerical 
simulations of Eq. (la) using the numerical schemes of Sec. 2. 

To study the shape of the stabilized soliton we define the error functional 

\\u(x,t)\ - \R p (x)\\\ 2 



\u(x,t)\\i 



(15) 
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Fig. 3. Same as Fig. 2 but with absorbing potential. 



and the error function E(t) as E(t) = min^ E^(t). Note that according to Eq. 
(5) the parameter /i parametrizes the shape of the Townes solitons in such a 
way that if \x increases the width of the soliton decreases and the amplitude 
increases. The error function E(t) measures the distance between the solution 
of Eq. (la) and the family of Townes solitons. We have computed numerically 
the function E(t) by evaluating the integrals in (15) for a large set of \x values 
and then choosing the minimum between all of them. 



In Fig. 5 we summarize our results when a Townes soliton is taken as initial 
data. We see, in Fig. 5(b) that E oscillates about a value of 0.02, with a 
maximum of about 0.03. This means that the solution of Eq. (la) can be 
approximated by a Townes soliton at all times with a maximum error in the 
norm of about 3% which is definitely small. In fact, this number could be 
even compatible with zero taking into account the number of approximations 
involved in our computation: the discretization of Eq. (la), the absorbing 
potential, the calculation of the integral (15) and the discrete set of fi values 
chosen to compute the minimum. 

When Gaussian initial data are taken (Fig. 6) we start with an optimal fitting 
whose error is about 14%. This error decreases with time to a value around of 
2% similar to the one in the Townes soliton case. This provides a quantitative 
support to the idea that Gaussian (or other) initial profiles evolve into a 
stabilized Townes soliton by ejecting a fraction of their amplitude in the form 
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Fig. 4. Results of numerical simulations of Eq. (la) showing stabilization of Gaus- 
sian initial data u(x,0) = (1/ ' ^pnw)e~ r l 2w with w = 1.09 for parameter values 
go = — 2-7T, gi = 8ir, £2 = 40. (a) Evolution of the width W(t) and of the amplitude 
max-j; The insets show pseudocolor plots of |tt(x,i)| 2 for different times, (b) 

Evolution of the norm II u II %. 



of radiation. 

Moreover, in Figs. 5(a) and 6(a) we see that the minimum points // m j n evolve 
according to what we expected. They oscillate in phase opposition respect to 
the width evolution of the solutions (see Figs. 3 and 4), because an increase 
of the width corresponds to a decrease of \i and viceversa (see Eq. (5)). 



3.2 Robustness of stabilized structures to parameter variations 



Finally, to have an idea of the range of parameter values for which stabilization 
is possible, we have made simulations when a Townes soliton is taken as initial 
data moving the g\ and Q parameters (#o — — 27r fixed). We have found that 
for Q G [20, 250] and g\ G [57r, 157r] the stabilization is possible although 
evolution of the width and of the amplitude shows different behaviors. In Fig. 
7 we plot the results for several parameter values. 
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Fig. 5. Analysis of the solutions of Eq. (la) for a Townes-soliton initial data, (a) 
Values of /x which minimize the error functional E„ as a function of time, (b) Error 
function E as a function of time. 




Fig. 6. Same as Fig. 5 but for Gaussian initial data. 



4 Conclusions 



In this paper we have developed an accurate numerical scheme to study the 
stabilization of solutions of the two-dimensional cubic nonlinear Schrodinger 
equation under variations of the nonlinear coefficient. We have also given nu- 
merical evidences in favor of the conjecture raised in Ref. [10] concerning the 
shape of stabilized solitons. Finally, we have shown that stabilization is possi- 
ble in a wide range of parameter values which is an interesting result for the 
possible applications of these structures to the fields of Bose-Einstein conden- 
sation and nonlinear Optics. 
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Fig. 7. Stabilization of u(x,Q) = R»,n(x) with fj, = 0.5 and g = -0.5 (W(0) = 1.09). 
Shown are the evolution of the width W(t) and of the amplitude max^ i)| for 
several parameter values, (a) go = —2tt, g\ = 6ir, Q = 20. (b) go = —2ir, g\ = 8ir, 
U = 250. (c) g = -2tt, gi = 6tt,Q = 80. (d) g = -2tt, 51 = 15tt, 17 = 150. 
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